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Abstract 

We present a model of mesoparticles, very much in the Dissipative Particle Dynamics 
spirit, in which a molecule is replaced by a particle with an internal thermodynamic degree of 
freedom (temperature or energy). The model is shown to give quantitavely accurate results 
for the simulation of shock waves in a crystalline polymer, and opens the way to a reduced 
model of detonation waves. 

1 Introduction 

Multimillion atom simulations are nowadays common in molecular dynamics (MD) studies. How- 
ever, the time and space scales numerically tractable are still far from being macroscopic, so that 
reduced models are of primary interest when multiscale phenomena are considered. In particular, 
the simulation of shock waves is a challenging task, involving very small time and space scales and 
large energies near the shock front, and much larger time and space scales and lower energies for 
the relaxation of the shocked materials, including the evolution of dislocations loops for example. 

Some reduced models for shock waves were proposed, for polycrystalin materials jH], crystalin 
materials with a projection of the dynamics onto a one-dimensional atom chain or resorting 
to mesoparticles with internal degrees of freedom pH] ■ This last approach seems to be the most 
promising and the most general one, and consists in replacing a complex molecule by a single 
particle. The introduction of an internal degree of freedom describing in a mean way the behavior of 
several degrees of freedom is reminiscent from Dissipative Particle Dynamics (DPD) models, which 
aim at describing complex fluids through some mesodynamics with some additional variables. 

DPD models, introduced in pHli have been put on firm thermodynamics ground in jg]. Some 
derivations from MD where proposed in a simplified case in 0, the more convincing general 
derivation being at the moment [J]. These studies motivate the modelling of the mean action 
of the projected degrees of freedom through some dissipative forces (depending on the relative 
velocities of the particles, so that the global momentum is conserved), balanced by some random 
forces. Ergodicity of the dynamics can be shown in some simplified cases Therefore, DPD 
dynamics are well established and motivated reduced models. 

Coarser models such as SPH (Smoothed particle hydrodynamics) [llLll^^j are routinely used to 
simulate shock waves at the hydrodynamic level, and can also be formulated in a DPD framework 
(the so-called Smoothed dissipative particle dynamics |H]). However, these models require the 
knowledge of an equation of state Eint = Eint{S,P) giving the internal energy as a function of 
entropy and pressure, for instance. Therefore, SPH-like models cannot be considered when the 
coarse-grained model is still at the microscopic level. 



It is also not possible to resort directly to the classical DPD models to simulate shock waves. 
Indeed, the dissipative and random forces arising in DPD are linked through some fluctuation- 
dissipation relation, using a local temperature. But when a shock wave passes, energy is transfered 
to the material, and the local temperature changes. Therefore, it is necessary to consider DPD 
models where the fluctuation-dissipation relation is not flxed a priori, but evolves depending on 
the physical events that have happened. DPD with conserved energy U are such models. 
Let us emphasize at this point that keeping thermal fluctuations in the microscopic models is of 
paramount importance to obtain the right relaxation proflles behind the shock front (HI El ■ 

We present here a dynamics strongly inspired by those models, and show that it provides an 
interesting mesoscopic model for the simulation of shock waves. To our knowledge, this is the 
flrst study on shock waves using stochastic dynamics of this kind. It also opens the way for an 
extension to detonation waves 52] , where exothermic chemical reactions are triggered as the shock 
passes, with the shock sustained and enhanced through the energy released. 

The paper is organized as follows. We flrst recall the usual DPD equations, and propose our 
mesoscopic dynamics. We then discuss the numerical implementation of the dynamics and present 
some simulation results. 



2 Dissipative Particle Dynamics with conserved energy 

All atom simulations are performed resorting to Newton's equations of motion. The correspond- 
ing microscopic systems are deterministic, Galilean invariant, and have some invariants, such as 
the total energy. While stochastic models are natural models to describe systems with reduced 
dynamics (since the information lost by the averaging process is modelled by some random pro- 
cess), it is however not clear that such a stochastic model can reproduce, even in a mean way, a 
deterministic dynamics with invariants. 

It turns out however that DPD models are stochastic dynamics which are Galilean invariant 
and preserve total momentum. Some reflnements were also proposed in order to conserve the total 
energy of the system, a model called 'DPD with conserved energy' (DPDE 

We consider a system of N particles in a space of dimension d, described by their positions 
(qi, . . . ,qN) and momenta {pi,...,pn), with associated mass matrix M = Diag(mi, . . . , tojv), 
interacting through a potential V. We assume for simplicity that the interactions between the 
particles are pairwise and depend only on the relative distances, so that V{q) = J2i<j 
Denoting by T the reference temperature and /? = 1/(^6?), the DPD equations read 

dqi — — dt 

mi 

dpt^ ^~W{rij)dt--/x'^ir,j){vi.j ■eij)eij + ^l^xin])dWije,j ^ ^ 



with 



7 > 0, r,j^\q,-q.j\ 



qj _ Pi _ Pj 



X a weight function (with support in [0,rc] where Tc is a cut-off radius), and where Wij are 
1-dimensional independent Wiener processes such that Wij = Wji. 

Notice that, since the dissipation term depends only on the relative velocities, the dynamics are 
globally GaHlean invariant. Besides, the total momentum is preserved. However, the total energy 
fluctuates, so that some reflnements in the model are required. Relying on the general DPD 
picture, DPD with conserved energy were introduced in pjl^. The idea is that the variations of 
the total mechanical energy 

H{q,p)^]^p^Mp + V{q) 

through the dissipative forces are compensated by some reservoir energy variable attached to each 
particle. Introducing an internal energy for each particle, the evolution of the internal energies 
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are constructed such that dH{q,p) + J2i — 0- associated entropy Si = s{ei) and an internal 
temperature can be also defined for each particle as 




For example, when the internal degrees of freedom are purely harmonic, T(e) = e/C„, where 
Cv is the specific heat at constant volume. More generally, this microscopic state law should be 
computed using MD or ab initio simulations. 



3 A simplified model for shock waves 

The model we consider is strongly inspired from DPD models with conserved energy U , so 
that all the properties of the usual DPD models with conserved energy can be straightforwardly 
transposed to this case. The derivation of the model is done as in QJIl]. 

The main differences here is that (i) we present the dynamics for particles of unequal masses, 
and (ii) do not project the dissipatives and random forces along the lines of center of the parti- 
cles. The generaHzation to particles of unequal masses is done by considering dissipation forces 
depending on the relative velocities, and not on the relative momenta. This is important if mix- 
tures composed of (say) two molecules are simulated, and each molecule is replaced by a single 
particle, whose mass is the total mass of the molecule. The dissipative and random forces could 
be projected as well to conserve angular momentum, but we restrict ourselves to the simpler and 
more general case when these forces are not projected, since we are only interested in Galilean 
invariance, and have in mind an extension to reduced models for reactive shock waves which 
do not necessarily preserve angular momentum, even if the dissipative and random forces are pro- 
jected. Such a model is also closer to the Langevin picture of the previous reduced models for 
shock waves O El ■ 

We finally neglect the thermal conduction here, since the contribution to the evolution of the 
internal energy arising from the dissipation forces is expected to be dominant in the nonequilibrium 
zone near the shock front. Heat diffusion plays a role only after the relaxation towards equilibrium 
in the shocked zone is achieved. 

The equations of motion for the system read: 



dqi = dt 
rrii 

dpi = ^ -VV'(ry) dt - JijX^{nj)vij dt + (Jijx{rij)dWij , 



(2) 



where x is still a weight function (with support in [0,rc] where is a cut-off radius), and Wij are 
now d-dimensional independents Wiener processes such that Wij = —Wji. The friction 7^ and 
the fiuctuation magnitude aij will be precised below. As for DPD models with conserved energy, 
the dynamics is postulated in a manner such that the total energy 

E{q,p,e) = H{q,p)+^e, 

i 

is preserved. The evolution of dH = — ^ ■ dei is inferred from ^ using Tto rule (see |1| for more 
details). Therefore, we consider the following dynamics: 

dqi = — dt 

dpi = ^ -W {n-j) dt - %jx'^{rij)vij dt + aijx{rij)dWij , 

jj^i (3) 

de^ = 2 I ^'^^'^vhv^l " ( ^ + ^ ) ^^^'''^^ I '^^ ^ "^'^ ^('^y)^ • dW^j, 

j, j^i \ \ I J / ) 
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with the following fluctuation-dissipation relation 13] 



It is then easily checked that measures of the form 

dpiq,P, e) = ie-^^(^'rt exp (^E ^ " Z^^^) ^^^^^^ (5) 

are invariant This measure expresses the fact that the translational degrees of freedom are 
distributed according to a classical Boltzmann statistics, whereas the internal energies are dis- 
tributed according to some free energy statistics. Since the total momentum Pq = J2iPi ^^'^ the 
total energy Eq = E{q,p, e) are also preserved by construction, the measure 

dpiq,p, e) = J—c-'^^ii'P) exp ( V 4^ - PeA Se^Eo Sp^p^ dqdpde (6) 

\^ Kb J 

is an invariant measure. 

If the dynamics is ergodic for the measure , it holds 



mt 



with Tkin = -fc^E^i 7^ = 7iJ2f=i ^) and (A) = J A{q,p) p{q,p,e) dqdpde. Notice that 
these relationships provide estimators for the local thermodynamic temperature f3~^/kB through 
the arithmetic average kinetic temperatures, and the harmonic average internal temperatures. Let 
us emphasize that a straightforward arithmetic average over the internal temperatures would give 
wrong results (the corresponding estimator being biased). The same result holds in the case when 
there are invariants of the dynamics (with the invariant measure lO) in the limit TV +00. 



4 A deterministic version of the model 

We intend here to introduce a deterministic version of our model, which allows to bridge the 
gap between a previous mesoscopic deterministic model JHl and the DPD framework for shock 
waves. The model proposed in JSI introduces damping forces on the position variable directly 
(and not on the momentum variables as would be expected) in order to preserve the Galilean 
invariance. Indeed, the damping terms in the momentum variable are considered to be of the 
form —j{vi — Vi), where Vi is a local average of the velocities around the particle, which makes 
the GaHlean invariance of the dissipated energy difHcult to preserve. If on the other hand the 
dissipation term in the momentum variable implies only pairwise velocity differences as for DPD 
models, the Galilean invariance follows immediately. The following equations of motion then mix 
the deterministic equations of motion of ^HI and the DPD philosophy: 

dqi 

dpi 

dei 

where T°j^^ is the average temperature in the kinetic degrees of freedom of particles i and j (for 
example, T°^' = {T°^'^ + T°^^)/2 with Tf^^ = 2p'l/kBdmi the kinetic temperature associated 



Pi 



dt 
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with particle i) and T'"* is the average internal temperatures of particles i and j (for example, 
Tj'"' = (T!/"* + Tj"')/2). The function uj is still a weighting function, and 7 determines the strength 
of the coupling. 

Notice that the dissipation term is in fact a dissipation term only when T°^* > T/"*, and an 
anti-dissipation term otherwise. This ensures that the internal and external (kinetic thus potential 
terms) energies equilibriate in all cases. However, the thermodynamic properties of such a model 
are less clear to state than for the previous stochastic model, and so, we stick to the model (PJ. 



5 Numerical discretization 



We use splitting formulas inspired from Recall that the integration of the equation of 

motion ^ is not straightforward since the dissipation terms depend on the relative velocities. 

We decompose (jSJ into elementary SDEs, and denote by (j)At the (stochastic) flow map for a 
time At. The elementary SDEs are the usual deterministic Newton part and the dissipation part, 
which read respectively 



dq - 
dp 



M-^pdt, 
-VV{q) dt 



and Vi < j, 



dpt = 
dp J = 

dei = 
dej = 



-li]X'^('rij)vij dt + ax{rij)dWij, 
-dpi, 

2 \2mi ^ 2m j J ' 
dCi. 



ri,3 

"diss, At 



{1 < i < i < N) the associated stochastic flow maps, an 

,1,2 ,N-1,N 



Denoting by 0Newton,At and 

approximation of (j>At of order 1 is 0At — ^dlss.At ° ■ ■ ■ ° <^di^,A" ° 0Newton,At- Notice that in 
practice it is difficult to resort to a second order accuracy sheme through a Strang splitting since 
this requires to keep track of the order the variables were updated in the first part of the time 
step, which is cumbersome when using Verlet lists. 

The Newton flow (/>Ncwton,At is approximated using a Velocity- Verlet scheme. For an approx- 
imation <&jfss At ^ of dissipation part, we flrst update the velocities at flxed internal 
temperatures using a Ver let-like algorithm as proposed in |16j : 



! 



pT'^' = p] + 7^l^,x\nM, - l<y^tx{nj) ur„ 



n+1/2 



n+l 
Pt =P: 



1 



1 



n+l 
ij 



-<T^Mx{n,)Ul 



where (?7i])i<i<j<Ar,n>o are independently and identically distributed standard gaussian random 
variables. Notice that the third and fourth updates can be rewritten in an expHcit form J^l- The 
energy is then updated as 



n+l 



1 {P^'f , 



2m, 



2m 1 



2m,- 




so that the total energy is indeed conserved by this step. Of course, this integration scheme could 
be refined, especially the dissipation part. We however tested several refinements and/or higher 
orders integrations in some model cases, but found no noticeable differences in the results. 



6 Application to shock waves 

Some numerical simulations of DPD models with conserved energy where proposed in |14[ [5] , 
but were concerned only with the computation of thermal conductivities. The corresponding 
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nonequilibrium states were stabilized using steady temperature gradients. The dissipation terms 
in the DPDE equations of motions were discarded, and only the diffusive part was retained. We 
present in this section profiles obtained from simulations of shock waves, for which the diffusive 
part of the dynamics can be discarded, but the dissipative part is of paramount importance to 
reproduce qualitative and quantitative features of all-atom shock waves. This situation is somehow 
complementary to the cases studied in ^^E], and, to our knowledge, was never considered before 
for some physical application. 

We consider the crystaUine polymer (PVDF) system of JHI, the corresponding reduced system 
being modeled by a two-dimensional (2D) triangular lattice of mesoparticles. Results for the 
all-atom model can also be found in |18) . 

The effective interaction potential between mesoparticles is a pairwise Rydberg potential of 
the form [iH] 



The parameters given by JHI were fitted to reproduce the stress in an uniaxial compression: 
A = 7.90, a = 0.185, ro = 5.07 A, e = 1.612 x lO^^o ^ = 64.03 x lO^^ kg/mol. We also choose 
a cut-off radius i?cut = 15 A for the pairwise interactions. The microscopic state law is obtained 
by assuming that is independent of the temperature: e = C„T, with here = 16 fee since 
we represent a three-dimensional molecule formed of 6 atoms by a 2D mesoparticle. In general, 
the heat capacity is a function of the temperature Ci, = Cu(T), and should be parametrized by 
equilibrium simulations. 

We use the simple weight function xif) = (1 ^ r/Rcut)^ if r > i?cut, x{i^) = otherwise, 
the cut-off radius i?cut being the same as the one used for the potential. Of course, many other 
weight functions could be used. We also set 7 = 1.5 x 10"^"' kg/s and At = 10~^^ s. In these 
preliminary tests of the model, the parameter 7 was varied to obtain a good agreement with the 
all-atom results. However, it is expected that 7 is Hnked to some physical quantity, such as the 
decay rate of the relative velocities autocorrelation in an all-atom simulation, and could therefore 
be estimated using some preliminary small equilibrium simulations. 

We first prepare an initial state according to the invariant measure JHJ. To this end, we 
sample independently the internal energies according to the measure exp(— /?e -I- s{e)/kB) = 
^-igC„/fcB exp(— /?e), and the initial configuration in phase-space by thermaHzing a lattice initially 
at rest, using a Langevin dynamics. In this study, the initial temperature is Tq = 300 K, and the 
edge of the triangles in the triangular lattice is a = 5.13 A. 

We then produce a shock using a piston at velocity Up = 3000 m/s. Figure presents the 
relaxation behind the shock front for the 2D triangular lattice of mesoparticles subjected to the 
dynamics l^l. The results are in good agreement with the all-atom results of pHl- In particular, the 
final temperature is very close to the all-atom value (whereas it is of course greatly overestimated 
by the mesoscopic dynamics without coupHng) , and the time required for the internal temperatures 
and kinetic temperatures to equilibriate is almost the time needed in all-atom studies. 

7 Conclusion and prospects 

We have proposed a simplified DPDE dynamics 0] , and shown its ability to simulate shock 
waves in a simple model case. This method allows a computational gain of one order of magnitude 
in space in the particular example considered, since 18 degrees of freedom are replaced by 2 degrees 
of freedom (assuming that the computational cost scales linearly with respect to the number of 
particles, which is the case here since we resort to short-ranged potentials). Besides, the time 
steps used can be taken larger as the time steps required by straightforward MD simulations since 
do not need to resolve the high-frequency atomic vibrations. 

Further quantitative agreement could be obtained by resorting to a less simplistic microscopic 
state law T = Tie), and parametrizing more systematically the dynamics. In particular, the 
friction coefficient 7 could be obtained from equilibrium all-atom simulations. Memory effects 
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Figure 1: (color online) Temporal evolution of the temperature of a thin slab of material as the 
shock runs trough it : mean kinetic temperature Tkin in the direction of the shock (intermediate 
curve, red), mean internal temperature Ti„t (lower curve, blue). The corresponding results when 
the coupHng with the internal degress of freedom is turned off are also shown (upper curve, black), 
and a cartoon representation of the all-atom result from JSl for the kinetic temperature Tkin is 
also plotted (dark dashed Hue). 
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could also be introduced by resorting to history-dependent friction terms (in analogy with the 
so-called generalized Langevin equations). 

More importantly, this study opens the way to a reduced model of detonation waves. For 
detonation waves, chemical reactions have to be treated explicitely and accurately in order to 
obtain the right velocities. Work is in progress to incorporate the corresponding chemical processes 
in the current DPDE framework |12j . 
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